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Abstract 

The Macaulay2 package PHCpack provides an interface to PHCpack, a general- 
purpose polynomial system solver that uses homotopy continuation. The main method 
is a numerical blackbox solver which is implemented for all Laurent systems. The 
package also provides a fast mixed volume computation, the ability to filter solutions, 
homotopy path tracking, and a numerical irreducible decomposition method. As the 
size of many problems in applied algebraic geometry often surpasses the capabilities 
of symbolic software, this package will be of interest to those working on problems 
involving large polynomial systems. 

1 Numerical homotopy continuation and PHCpack 

Many problems in applied algebraic geometry require solving, or counting the solutions of, 
a large polynomial or rational system. PHCpack is an interface to the program PHCPACK, 
one of several efficient polynomial system solvers that use numerical homotopy continuation 
methods [6]. 

The basic idea behind homotopy continuation is simple: to solve a polynomial system 
/(x) = 0, one first constructs a system g(x) = that is easy to solve and then constructs a 
homotopy, H(pc(t)) = (1 — t)#(x) + tf(x), in order to numerically track paths from known 
solutions of g (with t — 0) to the solutions of the target system / (with t = 1). 

Available since release 1.4 of Macaulay2 [3J, this package is motivated by [_5J and 
uses the data types defined by Leykin in NAGtypes .m2. The main function of the package 
allows a Macaulay2 user to solve a system numerically through a blackbox solver, where 
the creation of the start system and homotopy continuation is done behind the scenes. The 
package also provides a fast mixed volume computation and allows the user to filter solutions, 
to track solution paths explicitly, and to perform numerical irreducible decompositions. 
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In fact, the interface PHCpack offers access to most of the functionality of the software 
PHCpack, which has been serving as a development platform for many of the algorithms in 
numerical algebraic geometry [8J. Computations in this paper were done with phc ver- 
sion 2.3.61 (version 1.0 was archived in [9"). Since version 2.3.13, PHCpack contains 
Mixed Vol [2], and more recently added features are described in [10]. Note that PHCpack 
can solve Laurent systems, so the package includes a method to convert a rational system to a 
Laurent polynomial system. The underlying polyhedral methods perform well on benchmark 
problems; in many of those, the mixed volume is computed essentially instantaneously. 

Although PHCpack is open source, we follow the idea of OpenXM [7] and require only 
that the executable phc is available in the execution path of the computer. 

2 Numerical solutions of a polynomial system 

The main function, solveSystem, returns solutions of a system of polynomial or rational 
equations. Solutions are returned using data types from NAGtypes: a collection of Points 
which are approximations to all complex isolated solutions, or a WitnessSet for positive- 
dimensional components. The following system consists of 21 polynomial equations in 21 
unknowns, related to a Gaussian cycle conjecture [TJ §7.4, page 159] in algebraic statistics. 
The corresponding variety is zero-dimensional of degree 67. 

Macaulay2, version 1.4 

1 1 : CC [x_ 1 1 , x_ 12 , x_ 16 , x_22 , x_23 , x_33 , x_34 , x_44 , x_45 , x_55 , x_56 , x_66 , 

y_13 , y_14 , y_15 , y_24 , y_25 , y_26 , y_35 , y_36 , y_46] ; 

12 : system = {x_ll+2*x_12+2/3*x_16-l , x_23*y_13+5/2*x_12+ll/2*x_22 , 

x_33*y_13+x_34*y_14+ll/2*x_23, x_34*y_13+x_44*y_14+x_45*y_15, 
x_45*y_14+x_55*y_15+ (45/2) *x_56 , x_56*y_15+ (22/3) *x_16+ (45/2) *x_66 , 
82/7*x_12+17/2*x_22+12/5*x_23-l, x_34*y_24+( 14/11) *x_23+(12/5)*x_33, 
x_44*y_24+x_45*y_25+ (12/5) *x_34 , x_45*y_24+x_55*y_25+x_56*y_26 , 
x_56*y_25+x_66*y_26+ (82/7) *x_16 , (12/5) *x_23+ (282/5) *x_33+ (102/14) *x_34-l , 
x_45*y_35+ (282/5) *x_34+ (102/14) *x_44 , x_55*y_35+x_56*y_36+ ( 102/14) * x _45 , 
x_16*y_13+x_56*y_35+x_66*y_36, 10/l*x_34+205/16*x_44+(30/2)*x_45-l , 
x_56*y_46+ (205/16) *x_45+ (30/2) *x_55 , x_16*y_14+x_66*y_46+ (305/25) *x_56 , 
305/25*x_45+517/7*x_55+(89/3)*x_56-l, x_16*y_15+517/78*x_56+(89/3)*x_66, 
(450/21) *x_16+ (89/3) *x_56+ (293/19) *x_66-l}; 

13 : time solutions = solveSystem system ; 

using temporary files /tmp/M2-5331-2PHCinput and /tmp/M2-5331-2PHCoutput 
— used 0.056381 seconds 

14 : # solutions 
o4 = 67 

Solutions are returned as a list, each entry being of type Point, which includes diagnostic 
information such as the condition number and the value of the path-tracking variable t. This 
allows one to decide if a solution is "good" by using peek, suppressed here in the interest 
of space. Note also that the names of temporary input/output files allow the user to access 
the details of the entire phc computation, if desired. 

15 : solutions_0 
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o5 = {-3. 34446-1. 36293*ii, 2 . 1944+ . 682742*ii , -. 0664982- . 0038353*ii , 
-2 . 90447-1 . 02757*ii , - . 0074284+ . 306877*ii , . 018283- . 00390056*ii , 
-.0018297-. 0708941*ii, . 23468+ . 062941*ii , -. 13257- . 0064993*ii , 
. 00187754+ . 000036983*ii , . 083551+ . 0025807*ii , - . 00348342+ . 000364664*ii , 
12. 0202-34. 4693*ii, 14 . 2637- . 25899*ii , 6 . 77557+ . 029139*ii , 
5.38612-.346571*ii, 9 . 67526+ . 18591*ii , 8 . 34346- . 39709*ii , 
10. 7841-27. 2306*ii, 11 . 3312+ . 8239*ii , 20 . 0039+ . 37216*ii} 
o5 : Point 

The solutions can be further refined as necessary. To best illustrate refinement, consider 
the same system as above, but where the rational coefficients have been changed to larger 
rational numbers. Let newSolutions be the solutions of this modified system newSystem, 
which can be found in the Appendix (suppressed here in the interest of space). 

Solutions with coordinates below or above a given tolerance can be extracted by zeroFilter 
and nonZeroFilter, respectively. In the following example, we ask for solutions whose 12th 
coordinate is effectively zero (i.e., smaller than 10~ 19 ). Then, we confirm this by refining the 
answer to precision 64; notice that the 12th coordinate is now on the order of 1CT 67 . 

i9 : smallSolution = zeroFilter (newSolutions , 11 , 1 . Oe-19) 

09 = {{.0677823, -.386278, .0204925, -1.44743, .982877, -.366596, -.435274, 

.725281, -.422346, .0841728, .0218581, 2.23306e-20, 46.7882, -12.922, 
-70.411, 8.2731, -10.9958, 202.197, -43.8649, 306.199, 198.688}} 

110 : time smallerSolution=refineSolutions (newSystem, smallSolution, 64) 

— used 0.008859 seconds 

010 = {{.0677823, -.386278, .0204925, -1.44743, .982877, -.366596, -.435274, 

.725281, -.422346, .0841728, .0218581, -1.4308e-67, 46.7882, -12.922, 
-70.411, 8.2731, -10.9958, 202.197, -43.8649, 306.199, 198.688}} 

Note that when refining solutions, phc also recomputes input coefficients to a higher pre- 
cision, since rational coefficients may not always have an exact floating-point representation 
when the precision is limited. 

3 Mixed volume 

If the system has as many equations as unknowns, the mixed volume gives an upper bound 
on the number of isolated solutions with nonzero coordinates. For sufficiently generic coef- 
ficients, this bound is sharp. The function mixedVolume is illustrated below. 

111 : time mixedVolume (system ) 

using temporary files /tmp/M2-4281-6PHCinput and /tmp/M2-4281-6PHCoutput 

— used 0.011375 seconds 

011 = 75 

This polyhedral computation is faster than solving the system and provides an upper 
bound on the number of complex isolated roots in the torus. Computing the degree is much 
slower (and we note that it takes just as long to verify that the variety is zero-dimensional): 

112 : time degree ideal (system) 
— used 767.432 seconds 

012 = 67 
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While mixed volume counts solutions on the torus, one can also compute the stable mixed 
volume, which counts solutions with zero components as well, by using optional inputs to 
the method mixedVolume. phc offers additional functionality and flexibility, not all of which 
we can illustrate in this short note. Most interestingly, mixedVolume offers an option to 
use a start system, and creates a polyhedral homotopy from a random start system to the 
given system. The interested reader is referred to the documentation of the package for more 
information. 



4 Numerical irreducible decomposition 

Given a list of generators of an ideal /, the package can also compute a NumericalVariety 
with a WitnessSet for each irreducible component of V(I). The example below appears 
m[Tj. 

113 : CC[xll,x22,x21,xl2,x23,xl3,xl4,x24] ; 

114 : System={xll*x22-x21*xl2,xl2*x23-x22*xl3,xl3*x24-x23*xl4}; 

115 : V=numericalIrreducibleDecomposition(system) 
writing output to file /tmp/M2-5241-2PHCoutput 

calling phc -c < /tmp/M2-5241-3PHCbatch > /tmp/M2-5241-5PHCsession 

output of phc -c is in file /tmp/M2-5241-2PHCoutput 

. . . constructing witness sets . . . 

preparing input file to /tmp/M2-5241-7PHCinput 

preparing batch file to /tmp/M2-5241-9PHCbatch 

. . . calling monodromy breakup . . . 

session information of phc -f is in /tmp/M2-5241-10PHCsession 
output of phc -f is in file /tmp/M2-5241-8PHCoutput 
found 3 irreducible factors 

ol5 = A variety of dimension 5 with components in 

dim 5: (dim=5,deg=4) (dim=5,deg=2) (dim=5,deg=2) 
ol5 : NumericalVariety 

Witness sets are accessed by dimension: 

116 : WitSets=V#5; 

117 : w=first WitSets; 

118 : w#lslrreducible 
ol8 = true 

In the above example we found three components of dimension five. Let's verify the 
solutions. 

119 : R=QQ[xll ) x22 ) x21 ) xl2 ) x23,xl3,xl4 ) x24] ; 

120 : System={xll*x22-x21*xl2,xl2*x23-x22*xl3,xl3*x24-x23*xl4}; 

121 : PD=primaryDecomposition (ideal (system) ) ; 
i22: PD/print 

ideal (x23*xl4 - xl3*x24, x21*xl4 - xll*x24, x22*xl4 - xl2*x24, 

xl2*x23 - x22*xl3, xll*x23 - x21*xl3, xll*x22 - x21*xl2) 
ideal (xl3, x23, xll*x22 - x21*xl2) 
ideal (xl2, x22, x23*xl4 - x!3*x24) 
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As we see, the dimension and degree of each component agree with the numerical calcu- 
lation: 

123 : PD/dim 

023 = {5, 5, 5} 

124 : PD/degree 

024 = {4, 2, 2} 
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Appendix 

The polynomial system used in some of the examples in Section 2. It has the same support as the 
example system system, but the rational coefficients have been changed. It also has 67 solutions. 
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16 : newSystem = {(22531/300)*x_ll+(821/70)*x_12+(4507/210)*x_16-l, 

x_23*y_13+ (22531/300) *x_12+ (821/70) *x_22, 

x_33*y_13+x_34*y_14+ (821/70) *x_23 , x_34*y_13+x_44*y_14+x_45*y_15 , 

x_45*y_14+x_55*y_15+ (4507/210) *x_56 , 

x_56*y_15+ (22531/300) *x_16+(4507/210) *x_66 , 

(821/70) *x_12+ (140953/11025) *x_22+ (12325/504) *x_23-l, 

x _34*y_24+ ( 140953/ 1 1025) *x_23+ (12325/504) *x_33 , 

x_44*y_24+x_45*y_25+ (12325/504) *x_34 , x_45*y_24+x_55*y_25+x_56*y_26 , 
x_56*y_25+x_66*y_26+ (821/70) *x_16 , 

(12325/504) *x_23+ (282013/5184) *x_33+ (10231/1440) * x _34-l, 
x_45*y_35+ (282013/5184) *x_34+ (10231/1440) *x_44 , 

x_55*y_35+x_56*y_36+ (10231/1440) *x_45 , x_16*y_13+x_56*y_35+x_66*y_36 , 
(10231/1440) *x_34+ (205697/16200) *x_44+ (30529/2520) *x_45-l, 
x_56*y_46+ (205697/16200) *x_45+ (30529/2520) *x_55 , 
x_16*y_14+x_66*y_46+ (30529/2520) *x_56 , 

(30529/2520) *x_45+ (5175321/78400) *x_55+(897/35) *x_56-l , 
x_16*y_15+ (5175321/78400) *x_56+ (897/35) *x_66 , 
(4507/210) *x_16+ (897/35) *x_56+ (293581/19600) *x_66-l} ; 

17 : time newSolutions=phcSolve (newSystem) ; 

using temporary files /tmp/M2-3582-2PHCinput and /tmp/M2-3582-2PHCoutput 
— used 0.054471 seconds 

07 : List 

18 : # newSolutions 

08 = 67 
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